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Abstract 


This document is intended to detail the theoretical basis, equations, references and data 
that are necessary to enhance the functionality of commercially available Finite 
Element codes, with the objective of having functionality better suited for the 
aerospace industry in the area of composite structural analysis. The specific area of 
focus will be improvements to composite interlaminar fatigue and progressive 
interlaminar failure. Suggestions are biased towards codes that perform interlaminar 
Linear Elastic Fracture Mechanics (LEFM) using Virtual Crack Closure Technique 
(VCCT)-based algorithms [1,2]. All aspects of the science associated with composite 
interlaminar crack growth are not fully developed and the codes developed to predict 
this mode of failure must be programmed with sufficient flexibility to accommodate 
new functional relationships as the science matures. 


Nomenclature 


B 

Co 

C, 

Cu 

Cycle Jump 
da/dN 
(da/dN); 
(da/dN),1 
DCB 

ENF 


Number of blocks in block spectrum loading 

Paris Law constant, intercept 

Paris Law constant, stress ratio corrected; intercept. Mode I only 

Paris Law constant, stress ratio corrected; intercept. Mode II only 
Finite element solution load step from Pinay to Pmin but not a true load cycle 
Crack growth rate in length per cycle 

Crack growth rate in length per cycle, mode I only 

Crack growth rate in length per cycle, mode II only 

Double Cantilever Beam test for measuring mode I fracture toughness 
End Notch Flexure test for measuring mode II fracture toughness 
Energy Release Rate 

Fraction of max load Pmax to define the magnitude of each load block i 
Factor to account for R-curve effect where Gp = G), + f(a) 

Mode I static initiation fracture toughness 

Mode II static initiation fracture toughness 

Mode III static initiation fracture toughness 


Total critical energy release rate (fracture toughness) based on mixed mode 
law 


Mode I static initiation apparent fracture toughness with R-curve effect 


Mode I energy release rate at cyclic maximum load 


Gi —max 
Gin —max 
Gmax 
Gmin 


Gonset 


Rrocal 
Rappliea 


R-curve 


Mode II energy release rate at cyclic maximum load 

Mode III energy release rate at cyclic maximum load 

Total crack tip energy release rate at cyclic maximum load, Pmax 
Total crack tip energy release rate at cyclic minimum load, Prin 


Total crack tip energy release rate required to initiate growth from existing 
flaw 


Energy release rate threshold Paris Limit value 

Total crack tip energy release rate Gr = G, + Gy, + Gry 

Energy release rate threshold value 

Difference in maximum and minimum ERR AG = Gygax — Gmin 


Normalized cyclic energy release rate, Mode I only, where 
gi = G1—max/Gic 

Normalized cyclic energy release rate, Mode II only, where 
Gu = Gy-max/Grc 

Normalized cyclic energy release rate, Mode III only, where 
Gin = Gin-max/ Gir 


Stress intensity factor at maximum fatigue load 


Linear Elastic Fracture Mechanics 

Mixed Mode Bend test for measuring mixed-mode I & II fracture toughness 
Power law fit constant for Gons5¢; method 

Power law fit constant for Gons¢, method 

Fatigue load cycles 

Load cycles to accumulate first matrix damage 

Load cycles to initiation delamination from initial flaw 

Load cycles to grow delamination 

Total load cycles 

Progressive Damage Analysis 

Minimum applied load or displacement set in a fatigue load cycle 
Maximum applied load or displacement set in a fatigue load cycle 
Stress ratio R = Omin/Omax 

Local crack tip stress ratio estimated from Ry oc¢qi = «| Gromit) Oramay 
Externally applied stress ratio estimated from Ranpiied = Pmin/Pnax 


Changing interlaminar toughness with increasing crack length 


Pi 


Subcripts 


Cc 


Superscripts 


Virtual Crack Closure Technique 

Element length in crack growth direction 

Paris Law constant, slope 

Paris Law constant, slope. Mode I only 

Paris Law constant, slope. Mode II only 

Paris Law constant slope adjusted for stress ratio and mode mix 
Paris Law constant slope adjusted for stress ratio. Mode I only 
Paris Law constant slope adjusted for stress ratio. Mode II only 
Paris Law constant to correct intercept for stress ratio. Mode I only 
Paris Law constant to correct intercept for stress ratio. Mode II only 
Paris Law constant to correct intercept for stress ratio. Mode I only 
Paris Law constant to correct intercept for stress ratio. Mode II only 
Paris Law constant to correct slope for stress ratio. Mode I only 


Paris Law constant to correct slope for stress ratio. Mode II only 


Critical 
Mode I only 
Mode II only 
Mode III only 
element index 


Total 


increment number 


Paris Law constant, slope 


Introduction 


There exists a need in the aerospace industry to efficiently evaluate composite interlaminar 
crack growth under cyclic loading where loading consists of a complex spectrum of load 
components. This need is analogous to current analysis methods for metallic structure; however, 
interlaminar failure must deal with a three-dimensional delamination growth constrained between 
the plies and subjected to mixed mode interlaminar tension (mode I) and interlaminar shear (mode 
II) crack tip loading. Interlaminar crack growth is characterized using the DCB (Double Cantilever 
Beam) coupon for mode I; ENF (End Notch Flexure) test for mode II; and MMB (Mixed-Mode 
Bending) test for mixed-mode I/II cases (Refs.[3,4,5]). The Composite Materials Handbook-17 
(CMH-17) [6] defines three phases of fatigue delamination damage: “...the total cumulative life 
to failure, Nr, may be predicted by summing the lives for the onset of matrix cracking, Nw, 
delamination onset from this matrix crack, Np, and stable delamination growth to a finite 
acceptable size, Ng.” The total cumulative life to delamination failure is calculated in Eq. (1). For 
this article, “failure” is defined in the broadest terms to indicate a state of damage that exceeds a 
given design objective to be defined by the using organization. This document concerns only the 
onset and growth from a flaw, so the number of cycles associated with the matrix damage 
nucleation, Ny, is not addressed at this time. 


Nr = Ny + Np + Ne (1) 


Many thermoset composites tend to fail in interlaminar modes via brittle fracture. In comparison 
to ductile polymers, these materials are expected to be less sensitive to frequency or the path 
history between minimum and maximum load within the expected operating load spectrum. 
Consequently, fracture mechanics methods with linear scaling may be used to accommodate 
complex loading once the fatigue delamination propagation relationships (Paris Laws) have been 
characterized. The basic Paris Law is a power law function 


da 
dN 
where da/dN is the increase in delamination length per cycle and Ging is the maximum energy 


release rate at the crack front at peak loading. The factors Cp and exponents f£ were obtained by 
fitting the curve to the experimental data obtained from fracture tests. 


= C° GR as (2) 


There are three key modifications required for the basic interlaminar Paris Laws characterized 
by the DCB and ENF tests. These modifications account for 


(1) crack growth resistance mechanisms (such as fiber bridging, fiber delving, etc.) in the 
form of delamination growth resistance curves (R-curves) 


(2) stress ratio (e.g., R = Omin/Omax ) 
(3) mode mixity (G;;/Gr) 


Post-fatigue residual strength analysis requires interlaminar failure evaluation by robust 
progressive interlaminar failure algorithms. Code improvements in the following areas will greatly 
enhance the usefulness of these codes for aerospace structures. 


Specific Capabilities Needed 


Targeted improvements to commercially available interlaminar fatigue codes include: 


Cycle accumulation, Np, associated with damage onset as described in Ref. [7, 8] 
(Lower priority) 


Calculate Ng using efficient 3D progressive interlaminar fatigue growth algorithm based 
on VCCT 


o Include efficient Linear Elastic Fracture Mechanics (LEFM)-based algorithm 
similar to Ref. [9] 


o Use element post-critical load “ramping” to release constraints in a gradual 
manner. This is to account for a crack front traversing the mesh at an oblique 
angle [10] 


Enable various mixed mode I, II & II fatigue delamination growth laws 


o Interpolate along a constant da/dN contour based on one of several possible 
interpolation schemes such as Benzeggagh-Kenane (BK) Law [11], Reeder Law 
[12] or Power Law [12] 


o Interpolate along a constant G contour based on Ramkumar Law [13] 
co Interpolation based on laws proposed by Kenane [14] or Blanco [15] 


o Tabular inputs based on Mixed-Mode Bending (MMB) [16] fatigue Paris Laws 
and simple interpolation 


Implement Paris Law forms to include stress ratio influence, where R = 0 in/Omax 


Implement “Local Stress Ratio” acquired dynamically from crack tip energy release 
rates (ERR) 


o Calculate local stress ratio from crack tip ERRs, Rypocat = V Gr—min/Gr—max 


o Requirement to cycle to Pmin in “cycle jumps” to capture nonlinear effects or to 
update the local stress ratio after limited crack growth. Specifically, acquire p 


Local 


after a Pmin — Pmax cycle jump. Stress ratio is assumed to not vary significantly 
over short crack extension between cycle jumps. 


Damage accumulation under spectrum loading and load reversal 
Damage accumulation under intermixed fatigue and static crack growth 
Post-fatigue residual strength load cycle 


Post-processing of results and visualization of delamination-front contours 


The targeted improvements to the progressive VCCT static delamination include: 


Multi-element release with ramping 


Improved convergence algorithm 


Potential future enhancements that may be considered include: 
e Static and Fatigue Pristine Initiation (Calculation of Ny defined previously) 
e Static Crack Migration 
e Fatigue Crack Migration 


e Tri-linear strain softening law for simulating R-curve effects in static progressive 
analysis 


e Maintain compatibility for use for simulating in-plane damage and, specifically, 
interactions between in-plane cracks and delaminations 


The following three major sections address the theoretical background for suggested code 
improvements in the areas of fatigue delamination onset, fatigue delamination growth, and residual 
static (delamination) strength prediction. A fourth major section then provides a brief overview 
of potential future enhancements. 


Suggested Approach for Onset of Delamination Growth from a Singularity 
(discontinuity, matrix crack, existing delamination) 


This section addresses the cycles, Np, which are associated with the onset of delamination 
growth from a discontinuity, matrix crack, or existing delamination and are calculated following 
the procedures in Refs. [7,8 ]. The user should have the option to count cycles associated with 
damage onset, Np. However, the user should not be required to activate the algorithm. The damage 
onset cycles may be calculated for each constrained node pair along the crack front. Nodes along 
the initial crack front will have an accumulated total life associated with damage onset prior to 
entering into the crack growth phase. Damage onset calculations will be locked out once 
progressive fatigue delamination commences. A typical mode I R = 0.1 Gonser curve is shown in 
Figure 1. Such a curve can also be non-dimensionalized by Gic. A progressive Damage Analysis 
(PDA) code would need to accept mo and m; for each material, fracture mode, environment, and 
(potentially) a range of stress ratios R. One would anticipate that a curve similar to Figure 1 may 
be generated for mode II shear crack tip loading. The onset curves are fit as a power law to 
delamination onset data generated from a series of constant amplitude loads. The curves may be 
adjusted based on stress ratio or mode mix. Adjustments for mode mix and stress ratios for the 
progressive fatigue growth algorithm are discussed later in this document. Similar relationships 
may be derived for growth onset calculations. 
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Figure 1: Gonset curve from [7, 8] 


This section on delamination onset cycle counting is placed here based on the order of 
occurrence in an actual problem. However, the next section on progressive fatigue crack growth 
prediction, Nc, would be expected to be a more significant contributor to the damage life and 
would be a higher priority over counting cycles to delamination onset, Np. Additionally, 
establishing the progressive growth first, would aid in adding the onset cycles to the simulation in 
a second phase. 


Suggested Approach for Prediction of Interlaminar Fatigue Delamination 
Growth 


This section outlines suggested guidelines for progressive interlaminar fatigue crack growth 
life, Nc. The discussion is separated into subsections covering a 3D progressive fatigue algorithm, 
various forms of the Paris Law, Paris Law limits/thresholds, mixed-mode interpolation schemes, 
definition of a local/dynamic R-ratio, block spectrum algorithms, load reversal and negative stress 
ratio, and fatigue damage visualization. 


Suggested 3D Progressive Fatigue Algorithm 


The remainder of the specifications with respect to fatigue focuses on damage growth from the 
crack tip, based on LEFM and Paris Law data. One example is the VCCT-based progressive 
interlaminar fatigue analysis that was published in Ref. [9]. The analysis assumed crack growth 
in a brittle composite matrix calculated based on Paris law and stress ratio, R > 0. The algorithm 
took advantage of the brittle fracture characteristics of thermoset composites, and assumed that 
damage accumulation depends only on the maximum and minimum crack tip Energy Release 
Rates (ERRs). Crack growth was assumed to be independent of path between minimum and 
maximum load states and the case of load reversal was not considered. The simple algorithm of 
Figure 2 was implemented in a User Element (UEL) in the ABAQUS™ finite element code and 
was based on VCCT calculated energy release rates [10]. The “element” consisted of a center pair 
of nodes connected by a spring with a very high stiffness (essentially constrained). A series of 
unconstrained node pairs around the perimeter would sense an approaching crack front and would 
be used to calculate the crack opening displacement used in a VCCT calculation of energy release 
rate. Once the crack is calculated to grow across the element length, then the constrained node pair 
is released. 


VCCT calculations may be accomplished in other ways based on constrained node pairs, and 
the remainder of this document refers only to the “constrained nodes” with the understanding that 
other implementations may calculate energy release rates and release constraints in an analogous 
manner. In the algorithm of Ref. [9], the load step is ramped up to the applied maximum fatigue 
load and the load is held constant as the fatigue delamination is allowed to grow in the model, as 
shown in Figure 3. The crack tip energy release rate, Gina, and a crack growth rate, da/dN, are 
calculated for each node pair along the crack front based on a simple Paris Law, da/dN = Cy: 
any. Given the element length, Ax, and crack growth rate, the number of cycles required to 
separate a constrained node pair may be calculated, assuming self-similarity and constant mode 
mix during delamination extension. In one time-step increment, the algorithm will determine the 
minimum number of cycles, Nmin, to fail the most critical constrained node pair. The node pair 
with the fewest cycles to failure is released and all other nodes along the crack front accumulate 
damage based on the Njyin cycles in the current increment, AG cain: j= i oa wr (Njuin) : 
(da/dN); . The cycle is repeated in the next time increment, where damage accumulated in prior 
time increments is accounted for in the calculation of N,,;, for the current increment. At least one 
node pair is released per time increment while holding the load constant at Prax. 


G for all active 
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Figure 2: Paris Law linear fracture algorithm based on VCCT from [9]. 
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Figure 3: Ramp load to Pra, to acquire Gmax for Paris Law crack growth rate calculation [9]. 


The original code in Ref. [9] had several limitations. This algorithm was developed only for 2D 
crack growth and was not configured to account for mode mix. The code used an a priori stress 
ratio entered as an input value, and did not determine stress ratio from the crack tip energy release 
rates. Node pairs that had accumulated partial fatigue damage did not account for the damage as a 
change in spring stiffness until a sufficient number of cycles had been counted to fully release the 
constraining spring. 


Finite element code software vendors may offer a static VCCT interlaminar crack growth 
capability with both instantaneous release and “ramping” release after the crack tip reaches its 
critical energy release rate. Although the instantaneous release may be computationally efficient, 
a 3D crack front traversing a mesh at an angle will create a multitude of “corner nodes” as shown 
in Figure 4(A). Calculations at these corner nodes will result in a high energy release rate and the 
crack front may advance at artificially low load levels. Reference [17] shows how VCCT with 
instantaneous release may lead to artificially low predicted loads. However, the simulation may be 
improved if the “ramping” feature depicted in Figure 4(B) is implemented. During “ramping”, the 
post-critical load of the constrained node pair (or contact algorithm) follows an unloading curve 
encompassing the area determined by the energy release rate, G, as shown in Figure 4(C). This 
“ramping” is important for more precise crack growth predictions. Compared to the static crack 
growth prediction, those based on Paris Law have shown to be significantly more sensitive to 
variations in computed energy release rates. Therefore, a preferable implementation of a 3D 
fatigue delamination capability would incorporate a post-release unloading curve. Additionally, 
energy release rate calculations may be needed at intermediate crack tip locations between nodes. 
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Figure 4: Ramping to represent crack growth part way across elements, Ref. [18]. 


Suggested Forms of Paris Law to Include Stress Ratio Behavior 


The input syntax and crack growth algorithm must accommodate the basic interlaminar Paris 
Laws constants. Additionally, beyond the simple example shown in the previous sub-section, 
codes may accommodate parameters to account for R-curve behavior, stress ratio and mode mix. 
The input syntax must support various materials with Paris Laws fit to somewhat different equation 
forms. Table 1 contains three commonly used equation forms for crack growth rate. These equation 
forms are evaluated with the objective of deriving a unified crack growth law which contains 
sufficient flexibility to capture stress ratio effects for general loading and for materials yet to be 
characterized. Stress ratios between 0 and | are first considered, and a separate discussion is made 
later in this report concerning load reversal (R<0). 


Table 1: Paris Law Equation Forms 


Form Crack Growth Rate Equation Reference 
Walker Law da = =¢li- R) i K.. [19] 
AG da Gh 
ae 6 
a =|-c/ AG] [6] 
AVG da : 
i [A )-c{Ven—Jemn FT | po 


where R = 


Oomin _, |GT-min 
—min ~ {-f min 
Omax GT-max 


The Walker Law was derived for fatigue crack growth in aluminum panels. The equation 
accounts for the dependence of the crack growth rate, da/dN, on the stress ratio, R = Omin/Omax - 
Its development is detailed in Ref. [19]. The equation was originally presented in terms of 
maximum stress intensity factor, Kmax, and is converted to energy release rate here for the 
application to interlaminar fracture [21]. Subscripts “I” and “II” are added to account for mode I 
and II interlaminar crack growth. Mode III is assumed to be equal to mode H, for now, until 
significant data is available for mode III fatigue [12]. The mode I and II crack growth rates (Paris 
Law) are characterized with the DCB and ENF fatigue tests, respectively. Paris Law constants, C, 


and C,;, result from the conversion of Walker Law from the stress intensity factor to energy release 
rate. 


Cr = Gr Bu =Py Van =Vn (3) 


A convenient form of the Paris Law is found by normalizing the maximum energy release rate 


66 99 


by the initiation fracture toughness, G, or G,,. as discussed in Ref. [22]. The lower case “g” is 
introduced to represent the normalized form of the energy release rate, where g=G,,,, /G, and 
0<g<l. Hence, the Paris Law in its simplest form, with slope and intercept modified by stress 
ratio and R-curve effects, is given as Equation (4). 


Sor 
any (4) 
c= toeT] B' = B-(1-R) 


Units of C’ = length / cycle, 6’ is non-dimensional 


Accordingly, the Walker Law can be written in normalized form: 


=Cy° la > R)" ‘Sy la (5) 


The other two forms of Paris Law in Table 1 are converted to forms that include the maximum 
energy release rate, Ginax, during cyclic loading. 


AG form: (G) SOR ee (6) 
AVG form: (=) = C:[1—R) + Gnaxl? (7) 


The three Paris Law forms are captured by factors on Gmgy of (1—R)’, (1—R?), and 
(1 — R)*, corresponding to the Walker Law, the AG form and the AVG form, respectively. The 
following form of the Paris Law will capture all three: 


[S)=c-1-a-) -<f @ 


The three forms represented by Eq. (8) all assume the slope of the Paris Law is unchanged. This 
may be true for the class of thermoset composites where data exists in the literature. However, for 
other materials, the slope may change as a function of stress ratio. The following equation forms 
will add some flexibility in fitting data where the slope varies with stress ratio. 


(se acbey al o(Z)ech-er eh — w 


The document will proceed with the slope modified by 6’ = B - (1 — R)? and subscripts are 
added to represent mode I and II crack tip loading. 


(4) =C, aa y' ay 


da | an li (I-Ry" 
—]| =C,-\l-R”}’- 
(FF), chee y 


Implementation of the mixed mode functionality given in Equation (10) would provide a great 
early benefit to industry. 


(10) 


Forms of Paris Law to Include R-curve Behavior 


The application of a Paris Law normalized by static initiation fracture toughness may yield 
crack growth rates that are higher than crack growth observed in tests, since the normalized Paris 


law will not account for mechanisms that retard crack growth, such as fiber bridging, fiber delving, 
etc. R-curve effects can be eliminated from the Paris Law by either dividing DCB crack growth 
rate data by the static R-curve (resistance curve), or by limiting crack growth to lengths which are 
mostly void of R-curve influence. Equation (10) may be modified to account for delamination 
growth resistance mechanisms that will retard crack growth as a function of increasing crack length 
as discussed in Ref. [23]. R-curve effects are accounted for by normalizing G using the static 
resistance curve, Gir = Gjc * f(a). R-curve effects lead to an observed decrease in Paris law crack 
growth rates for typical co-cured interfaces. R-curve effects are dependent on the cyclic load levels 
under which these crack growth resistance effects are measured. Typically, R-curve effects are 
smaller under cyclic load levels compared to those measured under static crack growth conditions, 
as discussed in Refs. [24,25]. Consequently, normalization with the static R-curve, as in Equation 
(11), is approximate. Equation (11) is mathematically simple. However, programming a PDA code 
to spatially adjust the R-curve, f(a), as a function of lineal length from an arbitrary 3D crack front 
may be challenging. This functionality may be a point of future interest. 
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(11) 


Suggested Threshold and Paris Limit Inputs 


Fatigue crack growth rate curves are defined in three regions; region I is the regime where the 
crack growth rate transitions from the (notional) no-growth threshold, Gzu, to the linear Paris Law 
regime; region II is the linear Paris Law region and region III is the high growth rate region where 
the crack growth rate transitions from the upper Paris Limit, Gp;, to quasi-static crack growth. This 
section describes the expected code behavior approaching the threshold of growth, G7, and above 
the upper Paris Limit, Gp,. The Paris Limit, Gpz, marks the end of the linear Paris Law and is 
where crack growth transitions from subcritical crack growth under cyclic loading to quasi-static 
crack growth. The threshold of growth, Gry, marks the point at which the crack grows so slowly 
as to not be of practical engineering interest and resides below the end of the linear zone of the 
Paris Law. Two options are proposed for behavior in regions I and II, however an assessment of 
the accuracy of the choice of transition is not made. As before, the threshold and Paris limit are 
non-dimensionalized. Figure 5 defines Gry and Gp. 


Ge Ge. (12) 


0< 27 < Zp, <I 


Studies have shown the threshold of growth to be influenced by the stress ratio, as discussed in 
Ref. [20]. One practical approach to account for the threshold is to use Equation (4) and adjust the 
terms. For example, if the mode I threshold is measured in a DCB test at a stress ratio of R ~ 0, 
then an estimated threshold is calculated for R # 0 based on an equivalent crack growth rate. 
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Figure 5: Simple Paris Law with Gry and Gp; defined. 


The code may offer options on behavior in the transition regions. The behavior in the transition 
zones may follow no transition, a simple transition behavior or an advanced transition behavior, 
as defined in the following. 


Option 1: No Transition: This suggested approach ignores the threshold and Paris Limit 
behavior altogether. The crack growth rate follows the Paris Law equation except for the condition 
where the energy release rates reach the critical value for static crack growth. Many practical 
problems may not be greatly controlled by crack growth at the low and high fatigue crack growth 
regimes. No inputs are required for gry and gpr. 


Option 2: Simple transition: This suggested approach will calculate no growth for ERRs 
below gw and static release for ERRs greater than gp,. 


If g<g,,,, then AN=AN _... of lead node pair and Aa =0 


If g>g,,, then AN =1 and Aa = Ax (one element length). 


Option 3: Advanced transition: The Paris Law is modified with the factor from Ref. [23] 
[11 — (Gry /Gmax)?*)/(1 — (Gmax/Gc)”2)] where D; and D> are fit constants. See Figure 6. The 
non-dimensional form is [(1 — (gry/9max)°")/(1 — (Gmax/9c)”2)|. Option 3 is not typically 
employed in fatigue crack growth analysis and may be a future enhancement if there is a 
demonstrated need for accurate characterization in the transition regions. 
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Figure 6: Full-fatigue delamination characterization plot with region I and IT transitions [23]. 


Data from Ref. [7] suggests that thermoset composite laminates may not reach a true “no- 
growth” threshold and use of a Gry in analysis is provided as an engineering convenience for use 
with materials that appear to approach a threshold behavior. 


Suggested Mixed Mode I, II Paris Law Inputs with Interpolation Schemes 


Typically, interlaminar crack growth data is acquired under pure mode I and mode II loading 
using the DCB and ENF test methods, respectively. Optionally, MMB testing may be performed 
to evaluate crack growth under mixed mode loading at discrete mode mix conditions. This section 
considers algorithms to interpolate crack growth rates at mode mix intervals that lie between the 
mode mix ratios at which testing was completed. Following the de facto process for static crack 


growth, mode I and II crack growth is measured with the DCB and ENF test, respectively. The 
MMB test is used to determine a mixed mode law which is used to interpolate between modes I 
and II. Mode II toughness is often assumed to be equal to mode II toughness. Paris law fits for 
IM7/8552 tape (areal weight 190 gm/m’) are given in Figures 7a and 7b (See Refs. [16, 23]). 
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Figure 7: IM7/8552 Paris Laws with crack growth rate interpolation schemes. 


Energy release rates are given in absolute form (Figure 7a) and non-dimensional form (Figure 
7b). Interpolation may occur parallel to the “G”’ axis along a constant da/dN contour or parallel to 


the “da/dN” axis along a constant G contour. Interpolation may be between the pure mode Paris 
Laws curves or using intermediate Paris Laws curves from the MMB test. Figure 7a and 7b 
highlight the various options for interpolation. Crack growth thresholds and Paris Limits are 
included. 


Interpolation along constant G line using the Ramkumar Law [13] is given in Equation (14). 
Each mode J, II, II crack growth rate equation is of the Paris Law form given in Equation (4). 


de ; 
“= Cc, (g, y '+C'y (g, y "+C ny (1 y’ oe (14) 


dN 

This interpolation scheme assumes the accumulated damage in each mode is additive in the 
mixed mode condition. Although this equation is easy to use, it has not been supported by 
significant data. The two interpolation schemes along a constant da/dN contour, given as Eqs. (15) 
and (16), make the assumption that the mode mix law demonstrated under static growth translates 
to fatigue growth, thereby avoiding the requirement to perform separate MMB fatigue testing. 
These mode mix interpolation schemes have not been proven by significant data. A challenge 
with interpolation schemes along a constant da/dN is that da/dN may not be calculated explicitly. 
One approximate method is to forward calculate Gc~w for the target mode mix at two different 
da/dN values. Given two points, a new linear (approximate) mixed-mode Paris Law may be 
calculated. 


Interpolation along constant da/dN using Power Law: Solve for da/dN when gr= 1: 
G. 
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Following the two-point linear-fit process described above, both constant da/dN methods may 
result in a mixed-mode Paris Law. 


da 
dN 


Two additional mode mix laws that have been studied are discussed in Refs. [14, 15]. 
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From Kenane & Benzeggagh [14]: 
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Parameters C,, and 8. are additional fit parameters to be fit to MMB fatigue data. 


The advantage of equation (19) is that it does not assume monotonic behavior as a function of 
crack growth rate and may capture the unexpected mode mix behavior seen in Refs.15 and 16. 


Tabular Input: 


Another interpolation method requires measuring mixed-mode Paris Laws at various mode mix 
ratios [16] and performing a simple linear interpolation of C,,, and Bray: 


a = ‘(Gr anax 0 Su | _0,0,2, 0.4, 0.6,0.8, 1.0 (20) 
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Suggested Approach for Local Stress Ratio Acquired Dynamically by Simulation 


The energy release rate G computed at a local stress ratio, Riocai = V Gr—min/Gr—max always 
has to be compared to a Paris Law obtained for the same stress ratio in order to obtain the 
appropriate growth rate da/dN. The analysis load step, however, is dependent on the applied 
(external) stress ratio, Rupplied, Which may or may not be equal to the local stress ratio, Ryocai. 
Negative (local) stress ratios will be discussed in the next section. 


The most straightforward coding option for mixed mode Paris Law inputs is to not provide a 
stress ratio adjustment (i.¢., Rrocai= 0). In this scenario, the analyst will verify that the stress ratio 


for the mode I and II input Paris Laws is consistent with the simulation over the full crack growth 
range, or the user will make appropriate adjustments to the Paris Law constants external to the 
simulation (See Equation (4)). 


The second option is for the stress ratio to be calculated from the applied load and the 
assumption is made that the crack tip stress ratio trends with the applied stress ratio (Rrocai = 
Ru4pplied). The user will be required to confirm the compatibility of Paris Law constants in case the 
applied load results in a negative stress ratio. The code will require a method to identify the specific 
load component to define the stress ratio. 


The calculation of a local stress ratio was discussed previously (see Eq. 3) and should be made 
as the third option in the analysis setup. This feature may be important for problems where residual 
stress is important and the base (unloaded) state of the model is changing during the simulation 
(Riocai # RApplied). This option is likely a higher priority than option | or 2. In summary: 


e Option 1: No stress ratio adjustment or the stress ratio behavior is implied in the Paris 
Law constants. 


e Option 2: Stress ratio is calculated based on an external load component 


e Option 3: Stress ratio is calculated from the total ERRs 


Suggested Approach for Spectrum Effects 


Various methods are cited in the literature to account for load spectrum effects in metal fatigue. 
Cyclic loading damage phenomena such as crack tip blunting due to plasticity, hysteresis due to 
viscoelastic effects or other nonlinearities are assumed to not contribute to damage growth for 
interlaminar fatigue. Neglecting these nonlinear effects does not constitute a statement that such 
mechanisms are not active in all composite fatigue, only that the science and data is immature in 
this regard. Users of the code will need to assess the behavior of their material in question and 
ascertain if the LEFM methods proposed here are sufficiently applicable to their material. 


The influence of load spectrum on interlaminar crack growth is easily accounted for via the 
block spectrum loading damage accumulation algorithm cited in Ref. 26. This method is 
compatible with the energy release rate Paris Law methods that have become the norm for 
composite interlaminar crack growth. Simply stated, the total damage is the sum of the damage 
from multiple constant amplitude load blocks. No discussion is made here concerning the method 
by which the block spectrum is derived from the true spectrum. The crack is assumed to not grow 
an appreciable length through execution of the entire spectrum. The following discussion pertains 
to damage accumulation for interlaminar failure. The Paris Law given in the development is based 
on the mixed mode crack growth after adjustments have been made for stress ratio and R-curve. 
Each applied load will be in terms of maximum applied load, Pmax, and the corresponding 
minimum load, Pmin, for N cycles over B number of blocks. The model is initially loaded to Prin 
and then to Pa, to acquire crack tip energy release rates, Jmin and Gmax, respectively. 
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A convenient input syntax defines the individual blocks i as fractions fi 45, frin of the applied 


load, Prax [26]: 
eee Tak N, 
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Block spectrum damage accumulation methods assume that the material has negligible load- 
history effects. Load history effects would need to be characterized experimentally to evaluate the 
accuracy of the proposed approach. 


As the model is loaded from Prin to Pinax, and during every “cycle jump” from Prax to Pmin and 
back to Pmax, the energy release rate for each constrained node pair is acquired for every Dmin and 
Pinaxin the load cycle. 

i om fl ap 
- i tad is 
Pini = J iat? riage 
Ideally, the fatigue load step is configured such that min(p,in ) = Pmin and max(phax ) = 


Prax Where all other Divina Pinax in the block fall between Prin and Prax. The VCCT calculation 
provides the energy release rates: 
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(24) 
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For each block, the stress ratio, which is assumed to be constant during crack growth and 
updated at the next cycle jump, is calculated. As before, the user will have the option to calculate 
the stress ratio based on the applied load, R; = finin/fmax, or based on the local stress ratio (Eq. 


(24)). 


R. = Grin = Grin + Grin a9 Giron (25) 
i-Local Gi EG + Gime 
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The stress ratio and energy release rates are used to calculate the mixed mode Paris Law for 
each load block. Because the mode mix and stress ratio may change with each block, the mixed 
mode Paris Law constants, C an and B. are calculated uniquely for each block. The average 
crack growth rate is calculated and is subsequently used in the algorithm of Figure 2 to calculate 
the accumulated damage at each node pair along the crack front. 
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Suggested Approach for Load Reversals 


Metal fatigue is influenced greatly by the plastic zone around the crack tip and load reversal 
must be characterized by testing at a negative stress ratio. An analogous phenomena may occur in 
composite interlaminar fatigue. The analysis framework proposed here offers two mechanisms to 
account for load reversal: 


1. Provide Paris Law exponents compatible with the negative stress ratio to calculate the 
corresponding crack growth rate 


2. Treat the load reversal as a separate block loading and ignore loading history effects. In 
this scenario, load reversal at the crack tip will cause the shear (mode II and III) to 
reverse sign, however the mode I contribution will truncate at the zero load level as the 
sublaminate faces come in contact. 


A practical means to account for load reversal in a mixed mode problem is to divide the problem 
into a two-step block where tension will have mode I, IJ and II crack tip loading, and compression 
may have only modes II and III. Users will need to be aware that contact, preload or other 
nonlinearities may cause the zero energy release rate state to not coincide with the zero applied 
load. Currently, there is insufficient data for delamination growth under mode II conditions which 
allows a correlation of growth rates obtained from a fully reversed loading (R=-1) cycle with 
growth rates obtained from two cycles of a positive stress ratio (R>0). One may expect the 
mode I contribution to be the most significant part of the tension term of the growth rate expressed 
in Eq. (28). This is another area where the science may be immature and significant test-to-analysis 
correlation may be required to mature any code that accounts for stress ratio effects under reverse 


IN AVE IN Tension IN Compr ession 


Suggested Approach for Intermittent Static Growth, Pre-load and Residual 
Strength Load Cycle 


One would anticipate that a code could be configured to allow a progressive fatigue step to be 
interspersed with other static or fatigue load steps. Figure 8 shows one plausible loading scheme 
where a structure is preloaded statically (without damage growth), subjected to cyclic loading with 
fatigue damage growth, and subsequently evaluated in its post-fatigue damage state for residual 
static strength, at least for the delamination failure mode. The static preload will set the initial 
condition for the fatigue loading and the (minimum) crack tip energy release rates are recorded. 
Upon entering the fatigue step, the model is loaded to the maximum load. As describe previously, 


the maximum load may encompass a spectrum of loads, however the discussion will continue as 
if the structure is under constant amplitude cyclic loading. The discussion refers to singular applied 
loads, Pmin and Pmax, but these terms refer to a general set of loads and displacements applied to a 
structure at a “minimum” and “maximum” load condition. A stress ratio is calculated from the 
energy release rates at the maximum and minimum loads and is used in the crack growth rate 
calculations. The assumption is being made that the stress ratio does not change significantly over 
short lengths of growth. The maximum energy release rate and mode mix would be continually 
updated as the damage grows. In order to account for a changing stress ratio as damage grows, the 
minimum energy release rates are to be updated after each “cycle jump” back down to the 
minimum applied load. The cycle jump does not correspond to a physical load cycle, but is 
included to update the stress ratio to be used in the crack growth rate calculations. A simulation 
that is configured to run based on a predetermined stress ratio will not require a cycle jump and 
the user must confirm that the linear scaling assumptions are not violated by any non-linear 
behavior. In Figure 8, the minimum energy release rates Qmin(1), Zmin(2) ANd Ymini3), are different 
values which are used to update the stress ratios, Ry), Riz) and R,),. The number of cycles for 
every cycle jump may be predetermined or triggered based on some criteria. Progressive fatigue 
is expected to terminate after reaching the target life of the structure. 


A key question to answer is: Will the structure still meet its ultimate design loads after fatigue 
cycling for a full life of the structure (or inspection interval)? A final load step could be applied to 
simulate the target ultimate load with fatigue damage accumulated in the cyclic load step, e.g., to 
determine the residual static strength. An efficient progressive quasi-static interlaminar crack 
growth algorithm is important for accomplishing this analysis. Suggested capability is discussed 
in the following sections. 


Load 


© VCCT calculate gyi, 
@ VCCT calculate 1* g,,,, and R. Grow damage and count cycles Failure@ 
®@ Cycle Jump after Nc, cycles 

O Terminate progressive fatigue after target life, Nii. 


—— B max(i) 


Ria) 


——+ 8 max(i) 
Ri) | 


Damage | 
Growth | 


Damage 
Damage 
: *" Growth 


Growth | 


Piaii 
Time 
Incr. 
- >it Pie = > 
Static Progressive Residual 
Preload Fatigue Step Strength 
Step Step 


Figure 8: Possible loading scheme mixing fatigue and static damage growth. 


The simple algorithm suggested in Figure 8 implies that the most critical energy release rate 
coincides with the maximum load, Pmax. Likewise, 2min 1s assumed to coincide with Pin. A 
structure subjected to complex loading with contact and residual forces, however, may not have 
Lmax Coincident With Pnax and Ymin coincident with Pin. Additionally, the user may wish to analyze 
an array of loads for which the magnitude of certain loads is increasing while others are decreasing. 
Implementation of a complex algorithm to determine the applied load states corresponding to ginax 
and gmin may be too complex in the initial publication of a code capable of spectrum effects. A 
practical alternative is to interrogate values of Qmin and max on the load ramping from Piin to Pmax 
and report a warning to the user if the value of gma: does not coincide with Pmax or the value of gmin 
does not coincide with Prin. 


Suggested Approach for Post-Processing of Results, Visualization of 
Delamination Contours 


One would expect that a PDA code designed for interlaminar fatigue would have output 
consistent with the fatigue analysis. Data output may consist of contour plots to show: 
¢ Damage state as a function of cycles associated with a given increment 
* Current element state 
1) not active (not at crack front) 
2) opening but not released (e.g., blue color in shaded scale in a contour plot) 
3) releasing on ramp (e.g., green color in shaded scale in a contour plot) 
4) completely released (e.g., full red color in a contour plot) 
Other possible tabular output data may include the following: 
¢ Total cycles per increment 


* Damage area as a function of cycle count 


Suggested Improvements to Progressive Static Crack Growth 


The post-fatigue residual strength assessment is an important aspect in any structural fatigue 
evaluation. The predicted onset of damage using VCCT is a robust calculation. However 
progressive delamination growth predictions may be more challenging depending on the particular 
problem. This section identifies two key enhancement that will greatly improve the usefulness of 
currently available codes. 


¢ Adding multi-element release (within one time increment) at the iteration level with 
ramping as described in the previous sections. 


¢ Improvement of algorithm convergence. Specific approaches to this objective require 
further investigation; however, this topic remains a high priority. 


Specific Capabilities and Advanced Enhancements Suggested After Initial 
Tasks Above Have Been Implemented 


This section proposes advanced enhancements and capabilities that would improve the 
usefulness of VCCT-based interlaminar PDA codes. However, these enhancements are more 
appropriately implemented once the features described in prior sections have been implemented 
and fully vetted. The advanced topics that may be considered include: 


e Fatigue Pristine Initiation (Calculation of Ni defined in Ref. [7]) 
e Static Crack Migration analysis in Ref. [27] 
e Fatigue Crack Migration in Ref. [27] 


e Tri-linear strain softening law for simulating R-curve effects (progressive static 
analysis) in Ref. [17] 


e Maintain compatibility for use for simulating in-plane damage and specifically 
interactions between in-plane cracks and delaminations 
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